Leukemia Projection Tutorial

Here we will map an example dataset from three leukemias sequenced in our study. This entire tutorial should take less than 10 minutes to run

Setup

library(Seurat)
library(tidyverse)
library(symphony)
library(ggpubr)
library(patchwork)
library(RColorBrewer)

Warning: package ‘RColorBrewer’ was built under R version 4.1.2

Install package from github

## install dependencies that are not on CRAN
if(!require(BiocManager, quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("AUCell", "doMC", "BiocNeighbors"))
if(!require(devtools, quietly = TRUE)) install.packages("devtools")
devtools::install_github("jaredhuling/jcolors")

## install BoneMarrowMap package
devtools::install_github('andygxzeng/BoneMarrowMap', force = TRUE)

library(BoneMarrowMap)

Download reference object and UMAP model

Set projection folder, Download reference object and UMAP model

# Set directory to store projection reference files
projection_path = './'

# Download Bone Marrow Reference - 344 Mb
curl::curl_download('https://bonemarrowmap.s3.us-east-2.amazonaws.com/BoneMarrowMap_SymphonyReference.rds', 
                    destfile = paste0(projection_path, 'BoneMarrowMap_SymphonyReference.rds'))
# Download uwot model file - 221 Mb
curl::curl_download('https://bonemarrowmap.s3.us-east-2.amazonaws.com/BoneMarrowMap_uwot_model.uwot', 
                    destfile = paste0(projection_path, 'BoneMarrowMap_uwot_model.uwot'))

projection_path = './'
# Load Symphony reference
ref <- readRDS(paste0(projection_path, 'BoneMarrowMap_SymphonyReference.rds'))
# Set uwot path for UMAP projection
ref$save_uwot_path <- paste0(projection_path, 'BoneMarrowMap_uwot_model.uwot')

Visualize Bone Marrow Reference

If we want to visualize celltype labels or metadata from the BM Reference, we can create a Seurat Object from the symphony reference This will be memory efficient as it will not include gene expression counts, only the UMAP coordinates and the metadata including cell labels and sorting information

ReferenceSeuratObj <- create_ReferenceObject(ref)

DimPlot(ReferenceSeuratObj, reduction = 'umap', group.by = 'CellType_Annotation_formatted', 
        raster=FALSE, label=TRUE, label.size = 4)

We can also visualize broader cell type labels which may simplify interpretation and downstream analysis.

DimPlot(ReferenceSeuratObj, reduction = 'umap', group.by = 'CellType_Broad', 
        raster=FALSE, label=TRUE, label.size = 4)

We can visualize other annotations too, including cell cycle phase and lineage pseudotime estimates.

p1 <- DimPlot(ReferenceSeuratObj, reduction = 'umap', group.by = 'CyclePhase', raster=FALSE)
p2 <- FeaturePlot(ReferenceSeuratObj, reduction = 'umap', features = 'Pseudotime', raster=FALSE) 

p1 + p2

Load Leukemia scRNA-seq data.

We have confidently mapped leukemias spanning AML (incl. AMKL and AEL), B-ALL, MPAL, MDS, CML, MPN, and BDPCN, across sequencing technologies. However, we cannot be confident in mapping T-ALL due to a lack of thymus-specific reference data on T cell precursor stages. Further, our tool does not discriminate between normal vs malignant cells, which is an important consideration particularly within low-blast count chronic diseases.

As an example here, we are going to project scRNA-seq data from three diverse AML patients sequenced in our study.

• pt_17844: MLL-AF9 translocation with a purely mature leukemia cell hierarchy
• pt_17746: NPM1c + DNMT3A + TET2 with a GMP-dominant leukemia cell hierarchy
• pt_30886: Complex Karyotype with a primitive leukemia cell hierarchy + extensive erythroid involvement

For convenience, each patient sample was downsampled to 2500 cells.

curl::curl_download('https://bonemarrowmap.s3.us-east-2.amazonaws.com/ExampleQuery_AML_scRNAseq.rds', 
                    destfile = paste0(projection_path, 'ExampleQuery_AML_scRNAseq.rds'))

query <- readRDS('ExampleQuery_AML_scRNAseq.rds')
query

An object of class Seurat 
36601 features across 7500 samples within 1 assay 
Active assay: RNA (36601 features, 0 variable features)

Map the Query Data

Provide raw counts, metadata, and donor key. This should take <1 min Calculate mapping error and perform QC to remove low quality cells with high mapping error

# batch variable to correct in the query data, set as NULL if no batches in query
batchvar <- 'Patient'

# Map query dataset using Symphony (Kang et al 2021)
query <- map_Query(
    exp_query = query[['RNA']]@counts, 
    metadata_query = query@meta.data,
    ref_obj = ref,
    vars = batchvar
)

Normalizing
Scaling and synchronizing query gene expression
Found 2360 reference variable genes in query dataset
Project query cells using reference gene loadings
Clustering query cells to reference centroids
Correcting query batch effects
UMAP
All done!
Warning: Adding features not currently present in the object

In leukemia samples, the distribution of mapping error scores can vary broadly from sample to sample. In this context, we will want to threshold outliers with high mapping error on a per-sample basis. Typically, a threshold of 2, 2.5, or 3 MADs works well.

In some cases where sequencing depth is very low (e.g. older datasets from first-generation scRNA-seq protocols), a more stringent threshold of even 1.5 may be warranted to eliminate cells with low mapping quality

# Run QC based on mapping error score, flag cells with mapping error >= 2.5 MADs above median
query <- query %>% calculate_MappingError(., reference = ref, MAD_threshold = 2.5, 
                                          threshold_by_donor = TRUE, donor_key = batchvar) # threshold mapping error on a per-sample basis.

# Plot distribution by patient to ensure you are catching the tail
query@meta.data %>% 
  ggplot(aes(x = mapping_error_score, fill = mapping_error_QC)) + 
  geom_histogram(bins = 200) + facet_wrap(.~get(batchvar))

# Get QC Plots
QC_plots <- plot_MappingErrorQC(query)

# Plot together - If this is too crowded, can also just call "QC_plots" aloneto display one by one
patchwork::wrap_plots(QC_plots, ncol = 4, widths = c(0.8, 0.3, 0.8, 0.3))

This important step identifies a subset of cells with high mapping error from the query dataset that are either:

• not present within the reference, or
• have poor QC metrics (low RNA counts and low transcriptional diversity)

Sometimes, low quality cells may erroneously map to the orthochromatic erythroblast region as this cell type has very low transcriptional diversity. These low quality query cells do not have hemoglobin expression and are in fact mis-mapped; they will be flagged by the QC filter and excluded from cell type assignments.

Please adjust the MAD_threshold (typically between 1 and 3) based on the distribution of your dataset to identify the outliers with low quality and high mapping error scores. This will improve your classifications and any downstream composition analysis

# # Optional step - remove outliers with high mapping error
# query <- subset(query, mapping_error_QC == 'Pass')

Optionally, outlier cells with high mapping error can also be removed at this stage. For ease of integrating these mapped annotations with the rest of your analysis, we can choose to skip this step. If so, Final CellType and Pseudotime predictions will be assigned as NA for cells failing the mapping error QC threshold.

Cell Type Assignments

We will next use a KNN classifier to assign cell identity based on the 30 K-Nearest Neighbours from the reference map. Broader cell type labels will also be transferred automatically along with the precise cell type labels. This label transfer step should take <1 min for 10,000 cells and <5 min for 100,000 cells.

# Predict Hematopoietic Cell Types by KNN classification
query <- predict_CellTypes(
  query_obj = query, 
  ref_obj = ref, 
  initial_label = 'initial_CellType', # celltype assignments before filtering on mapping QC
  final_label = 'predicted_CellType'  # celltype assignments with map QC failing cells assigned as NA
) 

DimPlot(subset(query, mapping_error_QC == 'Pass'), reduction = 'umap_projected', group.by = c('predicted_CellType'), 
        raster=FALSE, label=TRUE, label.size = 4)

We can also visualize the broader cell type categories in case precise labels are too granular. These provide a simpler view of the data and can help guide cluster annotations if your dataset does not have many cells.

DimPlot(subset(query, mapping_error_QC == 'Pass'), reduction = 'umap_projected', group.by = c('predicted_CellType_Broad'), 
        raster=FALSE, label=TRUE, label.size = 4)

Pseudotime Annotations

We can also annotate each query cell based on their position along hematopoietic pseudotime. Query cells will be assigned a pseudotime score based on the 30 K-Nearest Neighbours from the reference map. Since our Pseudotime KNN assignments are performed in UMAP space (more accurate than KNN on harmony components), this step is very fast (< 10s)

# Predict Pseudotime values by KNN
query <- predict_Pseudotime(
  query_obj = query, 
  ref_obj = ref, 
  initial_label = 'initial_Pseudotime',  # pseudotime assignments before filtering on mapping QC
  final_label = 'predicted_Pseudotime'   # pseudotime assignments with map QC failing cells assigned as NA
)

# Visualize Hematopoietic Pseudotime in query data
FeaturePlot(subset(query, mapping_error_QC == 'Pass'), features = c('predicted_Pseudotime'), split.by = 'Patient') & 
  scale_color_gradientn(colors = rev(brewer.pal(11, 'RdBu')))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Visualize Projection Density

Now let’s visualize the density distribution of query cells across the hematopoietic hierarchy

# Set batch/condition to be visualized individually
batch_key <- 'Patient'

# returns a list of plots for each donor from a pre-specified batch variable
projection_plots <- plot_Projection_byDonor(
  query_obj = query, 
  batch_key = batch_key, 
  ref_obj = ref, 
  Hierarchy_only = FALSE, # Whether to exclude T/NK/Plasma/Stromal cells 
  downsample_reference = TRUE, 
  downsample_frac = 0.25,   # down-sample reference cells to 25%; reduces figure file size
  query_point_size = 0.2,   # adjust size of query cells based on # of cells
  saveplot = TRUE, 
  save_folder = 'projectionFigures/'
)

# show plots together with patchwork. Can also just call "projection_plots" object to display one-by-one
patchwork::wrap_plots(projection_plots, ncol = 3)

We can also set Hierarchy_only = TRUE to remove T/NK/Plasma/Stromal cells and focus solely on the hematopoietic hierarchy.

# Set batch/condition to be visualized individually
batch_key <- 'Patient'

# returns a list of plots for each donor from a pre-specified batch variable
projection_plots <- plot_Projection_byDonor(
  query_obj = query, 
  batch_key = batch_key, 
  ref_obj = ref, 
  Hierarchy_only = TRUE, # Whether to exclude T/NK/Plasma/Stromal cells 
  downsample_reference = TRUE, 
  downsample_frac = 0.25,   # down-sample reference cells to 25%; reduces figure file size
  query_point_size = 0.2,   # adjust size of query cells based on # of cells
  saveplot = TRUE, 
  save_folder = 'projectionFigures/'
)

# show plots together with patchwork. Can also just call "projection_plots" object to display one-by-one
patchwork::wrap_plots(projection_plots, ncol = 3)

Get Composition data for each donor

Here, to study the abundance of each cell type within each donor, I focus on cells that were classified with a KNN prob > 0.5 (that is, >50% of nearest neighbours from the reference map agree on the assigned cell type).

We can present this as a long table

query_composition <- get_Composition(
  query_obj = query, 
  donor_key = 'Patient', 
  celltype_label = 'predicted_CellType', 
  mapQC_col = 'mapping_error_QC', 
  knn_prob_cutoff = 0.5, 
  return_type = 'long')

query_composition 

Or as a wide table with counts of # of cells

query_composition <- get_Composition(
  query_obj = query, 
  donor_key = 'Patient', 
  celltype_label = 'predicted_CellType', 
  mapQC_col = 'mapping_error_QC', 
  knn_prob_cutoff = 0.5, 
  return_type = 'count')

query_composition 

Or as a wide table with proportion of each cell type within each donor

query_composition <- get_Composition(
  query_obj = query, 
  donor_key = 'Patient', 
  celltype_label = 'predicted_CellType', 
  mapQC_col = 'mapping_error_QC', 
  knn_prob_cutoff = 0.5, 
  return_type = 'proportion')

query_composition 

# Simple heatmap to visualize composition of projected samples
p <- query_composition %>% 
  # show celltypes present in >1% of total cells
  select(Patient, colnames(query_composition)[-1][colSums(query_composition[-1]) > 0.01]) %>% 
  # convert to matrix and display heatmap
  column_to_rownames('Patient') %>% data.matrix() %>% ComplexHeatmap::Heatmap()
p

Extra: Score AML cells for enrichment of LSC signatures

Here we will load genesets derived from functional studies of LSC+ and LSC- fractions in diverse AML patients (Ng Nature 2017), and LSC+ vs LSC- fractions in patients with MLL translocations (Somervaille Cell Stem Cell 2009). The former tends to correspond to more primitive LSCs while the latter corresponds to a distinct signature of mature promonocytic LSCs.

Further, we will load additional genesets derived from marker genes of leukemia cell populations defined in van Galen et al Cell 2019 and Zeng et al Nature Medicine 2022. Notably, the LSPC-Quiescent population is highly associated with functional LSCs and relapse by deconvolution (Zeng 2022).

curl::curl_download('https://raw.githubusercontent.com/andygxzeng/BoneMarrowMap_Extras/main/scAML_Differentiation_Stage_MarkerGenes.gmt', 
                    destfile = paste0(projection_path, 'scAML_Differentiation_Stage_MarkerGenes.gmt'))
AMLgenesets <- load_Genesets_gmt('scAML_Differentiation_Stage_MarkerGenes.gmt')

Read 20 items

AMLgenesets %>% summary()

                                   Length Class  Mode     
HSC_MPP_like__Top100_scAML         100    -none- character
LMPP_like__Top100_scAML            100    -none- character
Early_Lymphoid_like__Top100_scAML  100    -none- character
ProB_PreB_like__Top100_scAML       100    -none- character
MEP_MkP_like__Top100_scAML         100    -none- character
Early_Erythroid_like__Top100_scAML 100    -none- character
Late_Erythroid_like__Top100_scAML  100    -none- character
EoBasoMast_like__Top100_scAML      100    -none- character
GMP_like__Top100_scAML             100    -none- character
ProMono_like__Top100_scAML         100    -none- character
Monocyte_like__Top100_scAML        100    -none- character
cDC_like__Top100_scAML             100    -none- character
pDC_like__Top100_scAML             100    -none- character
LSC_vs_Blast_UP__Ng2016             47    -none- character
LSC_vs_Blast_DOWN__Ng2016           56    -none- character
MLL_LSC_UP__Somervaille2009        176    -none- character
MLL_LSC_DOWN__Somervaille2009      353    -none- character
LSPC_Quiescent__Zeng2022            61    -none- character
LSPC_Primed__Top100_Zeng2022       100    -none- character
LSPC_Cycle__Top100_Zeng2022        100    -none- character

# score genesets by AUCell and add to metadata
# typically I split the data into batches of 5k-10k cells to conserve memory
query <- score_Genesets_AUCell(query, genesets = AMLgenesets, nbatches = 1, ncores = 10, output = 'metadata')

[1] "Running AUCell scoring"
[1] "batch 1"

Warning: sparse->dense coercion: allocating vector of size 2.0 GiBQuantiles for the number of genes detected by cell: 
(Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).

    min      1%      5%     10%     50%    100% 
 507.00  613.99  740.00  838.00 1596.00 7032.00 

Using 10 cores.
Using 10 cores with doMC.
Genes in the gene sets NOT available in the dataset: 
    ProB_PreB_like__Top100_scAML:   2 (2% of 100)
    EoBasoMast_like__Top100_scAML:  1 (1% of 100)
    LSC_vs_Blast_UP__Ng2016:    1 (2% of 47)
    LSC_vs_Blast_DOWN__Ng2016:  6 (11% of 56)
    MLL_LSC_UP__Somervaille2009:    9 (6% of 160)
    MLL_LSC_DOWN__Somervaille2009:  6 (2% of 284)
    LSPC_Quiescent__Zeng2022:   6 (10% of 61)
    LSPC_Primed__Top100_Zeng2022:   9 (9% of 100)
    LSPC_Cycle__Top100_Zeng2022:    2 (2% of 100)
Garbage collection 179 = 115+25+39 (level 2) ... 
559.0 Mbytes of cons cells used (66%)
2131.0 Mbytes of vectors used (22%)

[1] "Finished Scoring"

query@meta.data

Using genes associated with functional LSC+ engraftment capacity (LSC_vs_Blast_UP__Ng2016_AUC) along with the LSPC-Quiescent signature help us to identify candidate LSCs within the data.

FeaturePlot(subset(query, mapping_error_QC == 'Pass'), 
            features = c('LSPC_Quiescent__Zeng2022_AUC', 'LSC_vs_Blast_UP__Ng2016_AUC'), 
            split.by = 'Patient', max.cutoff = 'q99', min.cutoff = 'q5') & 
  scale_color_gradientn(colors = rev(brewer.pal(11, 'RdBu')))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

It is no surprise that these are the most primitive cells earliest in Pseudotime. Note that this approach does not work for the MLL-AF9 samples wherein all leukemic cells are mature myeloid; in this context the MLL_LSC_Somervaille_2009_UP signature may help us identify functional “LSCs” capable of engraftment within the mature myeloid populations (typically late GMPs or Early ProMonocytes).

FeaturePlot(subset(query, mapping_error_QC == 'Pass'), 
            features = c('MLL_LSC_UP__Somervaille2009_AUC'), 
            split.by = 'Patient', max.cutoff = 'q99', min.cutoff = 'q5') & 
  scale_color_gradientn(colors = rev(brewer.pal(11, 'RdBu')))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

As expected, ProMonocytes within the MLL-AF9 sample exhibit enrichment for the MLL_LSC_Somervaille_2009_UP, in line with the mature myeloid origin of this patient’s disease.

Save projection results

This will save a csv file with the mapped annotations for each cell (mapping error scores, umap coordinates, predicted Cell Type, and predicted Pseudotime). We also have the option to save AUCell scores, provided that they are in the metadata.

# Save CellType Annotations and Projected UMAP coordinates
save_ProjectionResults(
  query_obj = query, 
  celltype_label = 'predicted_CellType', 
  celltype_KNNprob_label = 'predicted_CellType_prob', 
  pseudotime_label = 'predicted_Pseudotime', 
  save_AUCell_scores = TRUE,
  file_name = 'querydata_projected_labeled.csv')

Note on Downstream Analysis:

For downstream analysis, you can use the projected celltype labels to help annotate any leukemia cell clusters generated through unsupervised dimensionality reduction and clustering from individual patients. Sometimes, unsupervised analysis will yield clusters corresponding to different developmental states and other times it may yield clusters corresponding to distinct subclones within the patient. Along with tools like inferCNV for patients with known cytogenetic abnormalities, this can be integrated to visualize cellular hierarchies at the level of individual subclones.

Additional information from pseudotime projection and scoring of LSC-specific signatures can help identify candidate LSCs within the data.

For downstream analysis, let’s take the example of pt30886, with involvement of HSCs and multiple trajectories spanning Myeloid, Erythroid, and pDC differentiation:

pt30886 <- query %>%
  subset(., Patient == 'pt_30886') %>% 
  SCTransform() %>% 
  RunPCA() %>% 
  RunUMAP(reduction = 'pca', dims=1:20) %>%
  FindNeighbors(reduction = 'pca', dims = 1:20) %>% 
  FindClusters(resolution = 1, algorithm = 3)

Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 15690 by 2500
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 2500 cells

  |                                                                                                                                                                                                                          
  |                                                                                                                                                                                                                    |   0%
  |                                                                                                                                                                                                                          
  |=====================================================                                                                                                                                                               |  25%
  |                                                                                                                                                                                                                          
  |==========================================================================================================                                                                                                          |  50%
  |                                                                                                                                                                                                                          
  |===============================================================================================================================================================                                                     |  75%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================================================================================================| 100%

Found 84 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 15690 genes

  |                                                                                                                                                                                                                          
  |                                                                                                                                                                                                                    |   0%
  |                                                                                                                                                                                                                          
  |=======                                                                                                                                                                                                             |   3%
  |                                                                                                                                                                                                                          
  |=============                                                                                                                                                                                                       |   6%
  |                                                                                                                                                                                                                          
  |====================                                                                                                                                                                                                |   9%
  |                                                                                                                                                                                                                          
  |==========================                                                                                                                                                                                          |  12%
  |                                                                                                                                                                                                                          
  |=================================                                                                                                                                                                                   |  16%
  |                                                                                                                                                                                                                          
  |========================================                                                                                                                                                                            |  19%
  |                                                                                                                                                                                                                          
  |==============================================                                                                                                                                                                      |  22%
  |                                                                                                                                                                                                                          
  |=====================================================                                                                                                                                                               |  25%
  |                                                                                                                                                                                                                          
  |============================================================                                                                                                                                                        |  28%
  |                                                                                                                                                                                                                          
  |==================================================================                                                                                                                                                  |  31%
  |                                                                                                                                                                                                                          
  |=========================================================================                                                                                                                                           |  34%
  |                                                                                                                                                                                                                          
  |================================================================================                                                                                                                                    |  38%
  |                                                                                                                                                                                                                          
  |======================================================================================                                                                                                                              |  41%
  |                                                                                                                                                                                                                          
  |=============================================================================================                                                                                                                       |  44%
  |                                                                                                                                                                                                                          
  |===================================================================================================                                                                                                                 |  47%
  |                                                                                                                                                                                                                          
  |==========================================================================================================                                                                                                          |  50%
  |                                                                                                                                                                                                                          
  |=================================================================================================================                                                                                                   |  53%
  |                                                                                                                                                                                                                          
  |=======================================================================================================================                                                                                             |  56%
  |                                                                                                                                                                                                                          
  |==============================================================================================================================                                                                                      |  59%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================                                                                                |  62%
  |                                                                                                                                                                                                                          
  |===========================================================================================================================================                                                                         |  66%
  |                                                                                                                                                                                                                          
  |==================================================================================================================================================                                                                  |  69%
  |                                                                                                                                                                                                                          
  |========================================================================================================================================================                                                            |  72%
  |                                                                                                                                                                                                                          
  |===============================================================================================================================================================                                                     |  75%
  |                                                                                                                                                                                                                          
  |======================================================================================================================================================================                                              |  78%
  |                                                                                                                                                                                                                          
  |============================================================================================================================================================================                                        |  81%
  |                                                                                                                                                                                                                          
  |===================================================================================================================================================================================                                 |  84%
  |                                                                                                                                                                                                                          
  |==========================================================================================================================================================================================                          |  88%
  |                                                                                                                                                                                                                          
  |================================================================================================================================================================================================                    |  91%
  |                                                                                                                                                                                                                          
  |=======================================================================================================================================================================================================             |  94%
  |                                                                                                                                                                                                                          
  |=============================================================================================================================================================================================================       |  97%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================================================================================================| 100%

Computing corrected count matrix for 15690 genes

  |                                                                                                                                                                                                                          
  |                                                                                                                                                                                                                    |   0%
  |                                                                                                                                                                                                                          
  |=======                                                                                                                                                                                                             |   3%
  |                                                                                                                                                                                                                          
  |=============                                                                                                                                                                                                       |   6%
  |                                                                                                                                                                                                                          
  |====================                                                                                                                                                                                                |   9%
  |                                                                                                                                                                                                                          
  |==========================                                                                                                                                                                                          |  12%
  |                                                                                                                                                                                                                          
  |=================================                                                                                                                                                                                   |  16%
  |                                                                                                                                                                                                                          
  |========================================                                                                                                                                                                            |  19%
  |                                                                                                                                                                                                                          
  |==============================================                                                                                                                                                                      |  22%
  |                                                                                                                                                                                                                          
  |=====================================================                                                                                                                                                               |  25%
  |                                                                                                                                                                                                                          
  |============================================================                                                                                                                                                        |  28%
  |                                                                                                                                                                                                                          
  |==================================================================                                                                                                                                                  |  31%
  |                                                                                                                                                                                                                          
  |=========================================================================                                                                                                                                           |  34%
  |                                                                                                                                                                                                                          
  |================================================================================                                                                                                                                    |  38%
  |                                                                                                                                                                                                                          
  |======================================================================================                                                                                                                              |  41%
  |                                                                                                                                                                                                                          
  |=============================================================================================                                                                                                                       |  44%
  |                                                                                                                                                                                                                          
  |===================================================================================================                                                                                                                 |  47%
  |                                                                                                                                                                                                                          
  |==========================================================================================================                                                                                                          |  50%
  |                                                                                                                                                                                                                          
  |=================================================================================================================                                                                                                   |  53%
  |                                                                                                                                                                                                                          
  |=======================================================================================================================                                                                                             |  56%
  |                                                                                                                                                                                                                          
  |==============================================================================================================================                                                                                      |  59%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================                                                                                |  62%
  |                                                                                                                                                                                                                          
  |===========================================================================================================================================                                                                         |  66%
  |                                                                                                                                                                                                                          
  |==================================================================================================================================================                                                                  |  69%
  |                                                                                                                                                                                                                          
  |========================================================================================================================================================                                                            |  72%
  |                                                                                                                                                                                                                          
  |===============================================================================================================================================================                                                     |  75%
  |                                                                                                                                                                                                                          
  |======================================================================================================================================================================                                              |  78%
  |                                                                                                                                                                                                                          
  |============================================================================================================================================================================                                        |  81%
  |                                                                                                                                                                                                                          
  |===================================================================================================================================================================================                                 |  84%
  |                                                                                                                                                                                                                          
  |==========================================================================================================================================================================================                          |  88%
  |                                                                                                                                                                                                                          
  |================================================================================================================================================================================================                    |  91%
  |                                                                                                                                                                                                                          
  |=======================================================================================================================================================================================================             |  94%
  |                                                                                                                                                                                                                          
  |=============================================================================================================================================================================================================       |  97%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================================================================================================| 100%

Calculating gene attributes
Wall clock passed: Time difference of 26.35595 secs
Determine variable features
Place corrected count matrix in counts slot
Centering data matrix

  |                                                                                                                                                                                                                          
  |                                                                                                                                                                                                                    |   0%
  |                                                                                                                                                                                                                          
  |=====================================================                                                                                                                                                               |  25%
  |                                                                                                                                                                                                                          
  |==========================================================================================================                                                                                                          |  50%
  |                                                                                                                                                                                                                          
  |===============================================================================================================================================================                                                     |  75%
  |                                                                                                                                                                                                                          
  |====================================================================================================================================================================================================================| 100%
Set default assay to SCT
PC_ 1 
Positive:  HBB, AHSP, HBD, PRDX2, HBA1, CA1, HBA2, HBM, BLVRB, ALAS2 
       SLC25A37, HMBS, SLC4A1, IFI27, CA2, SLC25A39, HIST1H4C, MYL4, UROD, SLC2A1 
       HBG2, EPB42, RHAG, SMIM1, S100A6, DMTN, SELENBP1, UBAC1, HIST1H1E, BSG 
Negative:  CD74, MALAT1, HLA-B, HLA-DRA, HLA-DRB1, PMAIP1, BTG1, IFIT2, MX1, HLA-DPB1 
       HLA-A, AREG, SRGN, CXCR4, HLA-DPA1, CYBA, JCHAIN, FOS, SOX4, SAT1 
       EEF1A1, UBC, RPS27, OASL, HLA-C, RPL10, ZFP36L2, TMSB4X, HLA-DQB1, GPR183 
PC_ 2 
Positive:  GAPDH, SNHG29, RPLP0, C1QTNF4, PRSS57, GIHCG, PRDX1, RPS19, SERPINB1, AKR1A1 
       RPL13, MDK, GAS5, IFITM3, GSTP1, APOE, NAA38, ACTG1, VIM, TPI1 
       APOC1, CD34, BEX1, STMN1, NPDC1, SPINK2, CST3, FCER1A, SMIM24, RPL10 
Negative:  MALAT1, BTG1, CCL5, ISG20, MT-CO1, NKG7, HBM, HBA2, HBA1, CXCR4 
       HLA-B, SLC4A1, CA1, AHSP, MT-CO3, SLC25A37, CD7, TNFAIP3, CD3E, CD69 
       PRDX2, ADGRE5, CST7, XCL1, XCL2, NR4A2, UBC, RNF19A, SRGN, HBB 
PC_ 3 
Positive:  IFIT2, PMAIP1, CCL5, IFIT3, MT2A, OASL, CD7, TNFAIP3, MT-CO3, CCL3L1 
       CST7, ISG15, CD69, ISG20, CD247, IFIH1, CD3E, IFITM1, NKG7, MT-CO1 
       IFNB1, ZC3HAV1, XCL2, GZMM, XCL1, MT-CYB, IFIT1, HERC5, SPOCK2, TNF 
Negative:  CD74, JCHAIN, HLA-DRB1, HLA-DRA, HLA-DPB1, MX1, IGKC, PLD4, HLA-DPA1, LGMN 
       HLA-DQB1, SOX4, PLP2, AKAP13, CYBA, HERPUD1, CORO1A, AREG, LIME1, GPR183 
       IRF8, HLA-DRB5, HLA-DMA, SCT, SPIB, HLA-DQA1, PLAUR, SLC15A4, ITM2C, GPR65 
PC_ 4 
Positive:  HBB, S100A6, LYZ, ACTB, COTL1, CST3, S100A4, S100A10, HIST1H1E, TMSB10 
       PFN1, HBG2, IGKV3-15, IFI30, TMSB4X, TIMP1, FLNA, SH3BGRL3, EMP3, IGKV3-20 
       FTH1, VIM, ITGA2B, LGALS1, IGHV3-74, S100A9, TPI1, ANXA1, APOE, LMNA 
Negative:  C1QTNF4, IFIT2, PMAIP1, SPINK2, MDK, GYPC, SLC25A37, IFIT3, SLC4A1, IFNB1 
       GAS5, SMIM24, NPDC1, IGHM, CXCL10, HBA1, IFIH1, HBA2, SLC2A1, IFI27 
       HERC5, ZFP36L2, BEX1, SLC25A39, CRHBP, HEMGN, HMBS, ALAS2, EPB42, HBM 
PC_ 5 
Positive:  HBB, IFIT2, CCL3L1, CCL3, IFNB1, PMAIP1, IFIT3, CXCL10, HBG2, VIM 
       HIST1H1E, IGKV3-20, IGKV3-15, C1QTNF4, LYZ, IFIH1, CXCL11, IGHV3-74, CXCL8, HBM 
       NFKBIZ, S100A9, IFNL1, MPO, JUN, PTGS2, ATF3, BCL2A1, S100A8, CSF3R 
Negative:  APOC1, APOE, MYL4, EEF1A1, CD7, RPS27, CD3E, PVT1, FCER1A, IFITM1 
       NAA38, CNRIP1, BTG1, HLA-B, ATF7IP2, HES6, HERPUD1, ODC1, RNASE1, ATP5IF1 
       FAM178B, BLVRB, SDF2L1, CD247, RPS27L, TMSB4X, ITGA2B, MT-CO3, UROD, JCHAIN 
19:24:57 UMAP embedding parameters a = 0.9922 b = 1.112
19:24:57 Read 2500 rows and found 20 numeric columns
19:24:57 Using Annoy for neighbor search, n_neighbors = 30
19:24:57 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
19:24:57 Writing NN index file to temp file /var/folders/2b/f9nktt6d0wzfcpxv5xxm2n900000gp/T//RtmpGCyUkO/file11ec27047390
19:24:57 Searching Annoy index using 1 thread, search_k = 3000
19:24:57 Annoy recall = 100%
19:24:58 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
19:25:00 Initializing from normalized Laplacian + noise (using irlba)
19:25:00 Commencing optimization for 500 epochs, with 103452 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
19:25:03 Optimization finished
Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2500
Number of edges: 91096

Running smart local moving algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8258
Number of communities: 16
Elapsed time: 0 seconds

DimPlot(pt30886, reduction = 'umap', group.by=c('seurat_clusters', 'predicted_CellType_Broad'), label = T)

FeaturePlot(pt30886, 
            features = c('LSPC_Quiescent__Zeng2022_AUC', 'LSC_vs_Blast_UP__Ng2016_AUC', 'HSC_MPP_like__Top100_scAML_AUC'), 
            max.cutoff = 'q99', min.cutoff = 'q5', ncol = 3) & 
  scale_color_gradientn(colors = rev(brewer.pal(11, 'RdBu')))

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Based on this we can see multiple transcriptional clusters within the data, including some non-leukemic clusters of T cells, B cells, and Plasma cells. Cluster 3 contains HSC/MPP-like cells most enriched for LSC signatures.

Let’s see what each cluster looks like when projected onto BoneMarrowMap:

DimPlot(pt30886, reduction = 'umap_projected', group.by=c('seurat_clusters', 'predicted_CellType_Broad'), label = T)

pt30886@meta.data %>% 
  # compare clusters against broad cell types assignments 
  group_by(predicted_CellType_Broad, seurat_clusters) %>% 
  summarise(count = n()) %>% ungroup() %>% 
  # get proportions of cells within each cluster assigned as a certain celltype 
  group_by(seurat_clusters) %>% mutate(prop = count / sum(count)) %>% 
  pivot_wider(id_cols = 'predicted_CellType_Broad', names_from = 'seurat_clusters', values_from = 'prop') %>% 
  # fill in zeros and format for heatmap 
  replace(is.na(.), 0) %>% column_to_rownames('predicted_CellType_Broad') %>% 
  # create a basic heatmap to visualize cell types vs clusters
  pheatmap::pheatmap()

`summarise()` has grouped output by 'predicted_CellType_Broad'. You can override using the `.groups` argument.

This illustrates how annotations from BoneMarrowMap can allow us to analyze clusters of leukemic cells identified from unsupervised analysis. As shown in the heatmap above, each cluster from this patient actually corresponds to a distinct state along the hematopoietic hierarchy, with multiple clusters corresponding to erythroid progenitors. Overlaying copy number alteration calls from inferCNV or expressed mutations can help us identify the differentiation patterns of distinct subclones within this patient.

Finally, for cohorts with many patients composition analysis can be performed to understand how leukemia cell hierarchies vary with patient characteristics and therapy response / relapse. I hope this tutorial was useful and feel free to comment with any questions you may have around projection of leukemia cells or downstream analysis.